#include 'model_R2.h'
model_R2::model_R2(int iMode,double *H_profile, int L)
// Construct model_R2 from a height profile
{
	int i;
	Mode = iMode;
	Alpha = new double[Mode];
	Beta  = new double[Mode];

	H_profile = NULL;
	if H_profile == NULL
	{
		for(i=0;i<Mode;i++)
		{
			Alpha[i]  = 0;
			Alpha2[i] = 0;
			Beta[i]   = 0;
			Beta2[i]  = 0;
		}
	}
}


model_R2::~model_R2()
{
	delete [] Alpha2;
	delete [] Beta2;
}

double model_R2::set_state(double *iAlpha, double *iBeta)
{
	for(int i = 0;i<Mode;i++)
	{
		Alpha2[i] = iAlpha[i]*iAlpha[i];
		Beta2[i]  = i
double model_R2::getVarR2()
{
	return VarR2;
}

double model_R2::getMeanR2()
{
	return MeanR2;
}

model_R2::update(double w, double dt)
{
	int j;
	double temp1,temp2,temp3,temp4,temp5;
	double VarAlpha2,VarBeta2;

	nu = interp1(TAB_R_DEP_RATE,TAB_R_NU,w);
	sigma2 = interp1(TAB_R_REP_RATE,TAB_R_SIGMA2,w);
	
	// Update mean(R^2)
	MeanR2 = 0;
	for(j=1;j<=Mode;j++)
	{
		temp1 = sigma2/(2*nu*j*j);
		temp2 = exp(-2*nu*j*j*dt);
		Alpha2[j] = temp+(Alpha2[j]-temp)*temp2;
		Beta2[j]  = temp+(Beta2[j]-temp)*temp2;
		MeanR2 += (Alpha2[j]+Beta2[j])/(2*PI);
	}
	// Update var(R^2)
	VarR2 = 0;
	for(j=1;j<=Mode;j++)
	{
		temp1 = exp(-2*nu*j*j*dt);
		temp2 = exp(-4*nu*j*j*dt);
		temp3 = sigma2/(nu*j*j);
		temp4 = temp3*(2*temp2-2*temp1);
		temp5 = 0.5*temp3*temp3*(temp2-2*temp1+1);
		VarAlpha2 = -Alpha2[j]*temp4+temp5;
		VarBeta2  = -Beta2[j]*temp4+temp5;
		VarR2 += VarAlpha2+VarBeta2;
	}
}

